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Chapter 1 
Tests of Classical Gravity with Radio Pulsars 


Zexin Hu, Xueli Miao, Lijing Shao 


Abstract Tests of gravity are important to the development of our understanding of 
gravitation and spacetime. Binary pulsars provide a superb playground for testing 
gravity theories. In this chapter we pedagogically review the basics behind pulsar 
observations and pulsar timing. We illustrate various recent strong-field tests of the 
general relativity (GR) from the Hulse-Taylor pulsar PSR B1913+16, the double 
pulsar PSR J0737—3039, and the triple pulsar PSR J0337+1715. We also overview 
the inner structure of neutron stars (NSs) that may influence some gravity tests, 
and have used the scalar-tensor gravity and massive gravity theories as examples to 
demonstrate the usefulness of pulsar timing in constraining specific modified gravity 
theories. Outlooks to new radio telescopes for pulsar timing and synergies with other 
strong-field gravity tests are also presented. 


1.1 Introduction 


As one of the four fundamental interactions, gravity is the most mysterious to study. 
From Newton’s law of universal gravitation to Einstein’s general relativity (GR), in 
the development of gravity theory, experiments and observations have been playing 
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an important role in overturning or verifying different gravity theories. Currently, 
GR has passed all the tests with flying colors [118], and it seems that GR is the 
tour de force among all gravity theories. However, in GR, there is a fundamental 
difficulty in quantization as well as in describing spacetime singularities, suggesting 
that GR could not be a complete theory. It might be an effective field theory below 
some energy scale [113]. So, there is strong motivation to find deviations from GR 
predictions and to search for other possible alternative gravity theories. People have 
tried to construct quantum gravity theories to solve the quantization problem in GR, 
such as the string theory [40]. People have also used alternative gravity theories to 
explain the missing mass in galaxies and the accelerating expansion of our Universe 
without introducing the mysterious concepts of dark matter and dark energy [57, 
19]. To this end, we need to use all kinds of observations to probe different aspects 
of gravity and use measurement results to test gravity theories. 

People usually classify gravity tests based on the studied systems in different 
gravity regimes [117, 3, 58]. In astrophysics, for systems of astronomical scales, 
Wex [117] introduced the following four gravity regimes: (i) quasi-stationary weak- 
field regime, (ii) quasi-stationary strong-field regime, (iii) highly-dynamical strong- 
field regime, and (iv) radiation regime. In a quasi-stationary weak-field regime, the 
velocity of masses is much less than the speed of light c, and spacetime is close 
to a Minkowski spacetime. Our Solar system belongs to this regime. In a quasi- 
stationary strong-field regime, the velocity of masses is also much less than c, but 
one or more masses have strong self-gravity, which causes the neighborhood space- 
time to deviate significantly from a Minkowski one. A well-separated binary system 
containing one or two compact bodies belongs to this regime. In a highly-dynamical 
strong-field regime, the velocity of masses is close to c and spacetime is strongly 
curved. The merging of compact binaries reflects the gravity of such a regime. A ra- 
diation regime is the propagation regime of gravitational waves (GWs). Gravity tests 
from the Solar system have provided strict restrictions on the quasi-stationary weak- 
field regime [118], while the GW detection from LIGO/Virgo/KAGRA probes the 
radiation regime and provides opportunities to analyze the highly-dynamical strong- 
field regime [93]. To study the behavior of gravity in a quasi-stationary strong-field 
regime, we need to have one or more compact bodies in a well-separated system. 
Binary pulsar systems with good timing precision provide precise measurements of 
gravity in such a regime. 

In this chapter, we will provide a pedagogical introduction to pulsar timing, as 
well as some aspects and new developments of using it in testing classic (namely, 
not quantum) gravity theories. In the next section, we overview the basics of pulsar 
observations and the theoretical background for pulsar timing. We use the famous 
Hulse-Taylor pulsar PSR B1913+16, the double pulsar PSR JO737—3039, and the 
triple pulsar PSR JO337+1715 as primary examples for testing GR in Sec. 1.3. Af- 
ter that, we discuss the inner structure of neutron stars (NSs) in Sec. 1.4 which in 
some theories will affect the binary orbital dynamics. A couple of representative 
alternative gravity theories are tested against pulsar timing data in Sec. 1.5. Finally, 
Sec. 1.6 summaries the chapter. More reviews on the topic of using radio pulsars for 
gravity tests can be found in Refs. [106, 117, 78, 100, 97, 58]. 
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1.2 Pulsar Observation and Pulsar Timing 


In 1967, the first pulsar PSR B1919+21 was discovered [50]. The prevailing view of 
pulsars is that they are stably spinning NSs with a high-intensity magnetic field [45]. 
In a pulsar’s open magnetic field line region, the charged particles moving along 
magnetic field lines can produce curvature radiation, creating a radiation cone cen- 
tered on the magnetic axis [73]. Generally, a pulsar’s magnetic axis and spin axis 
do not coincide, so if the radiation region of a NS sweeps past Earth’s radio tele- 
scopes when the NS rotates, one can detect a series of periodic pulse signals on the 
Earth. In this setting, a rotating NS seems like a lighthouse, that is why we also call 
it the “lighthouse model.” Pulsars, as rapidly rotating NSs, provide an ideal labora- 
tory in astrophysics for fundamental physics. The masses of NSs are typically larger 
than that of our Sun, but their radii are only about 10 km. Therefore NSs represent 
an extremely dense stellar environment. Inside NSs there exists high matter density, 
which can reach and exceed nuclear density, ~ 10!> g/cm, as well as high pressure, 
~ 10°%erg/cm?. NSs also have high magnetic field ranging from 108G to 10°G. 

Benefiting from the development of telescopes, we have detected more than 3000 
pulsars up to now. Figure |.1 shows the population of pulsars by the distribution of 
spin period (P) and spin period derivative (P) [79]. In the figure the red circle marks 
pulsars in binary systems, which make up about 10% of detected pulsars. These 
pulsars’ spin periods range from 1.4ms to 24s. According to the distribution of P 
and P, one can classify pulsars into two groups: normal pulsars which correspond to 
the distribution in the upper right part in Fig. 1.1, and millisecond pulsars (MSPs) 
which correspond to the distribution in the bottom left part. For MSPs, as the name 
suggests, the spin periods are mostly less than 10 ms and the spin period derivatives 
range from 10~!8ss~! to 10-72 ss~!. The small P of MSPs means that they have 
more stable spin periods than normal pulsars on long timescales, and the stability 
of some MSPs can revive the precision of atomic clocks in long timescales [72, 
79, 53]. The fast rotation of MSPs is caused by the “recycling process” where the 
pulsars can spin up by a stable mass transfer from the companion stars [108]. So, 
MSPs are mostly found in binary systems, as clearly seen in Fig 1.1. Fully recycled 
pulsars would have millisecond spin periods and locate in very circular orbits. Some 
pulsars could undergo mild recycling processes and only spin up to = 10 ms, and 
their orbits would be eccentric. The stability of spin periods of MSPs makes the 
measurements of them highly precise. As discussed earlier, the motion of a MSP 
in a binary system probes the quasi-stationary strong-field regime of gravity. So 
detecting a stably rotating MSP in a binary system allows us to test gravity in a 
quasi-stationary strong-field regime. 

The radiation of a pulsar can be detected in multiple frequency bands. For exam- 
ple, the Crab pulsar, a young pulsar discovered in the Crab nebula, radiates pulses 
from the radio bands to y-ray bands [105, 75]. The discovery of the Crab pulsar in 
the Crab nebula also demonstrates that a rotating NS is a product of a supernova. 
While other bands also provide important information, study of pulsars is mainly 
benefited from the radio band. Pulsars usually have a power-law radio spectrum and 
the spectrum index is normally from —2 to —1.5 [73]. 
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Fig. 1.1: The period—period derivative diagram for pulsars [79]. The black dots 
are pulsars with pulsed emission in the radio band and the red circles mark binary 
systems. We also denote “AXP” for systems that are anomalous X-ray pulsars or 
soft y-ray repeaters, “RRAT” for systems that are pulsars with intermittently pulsed 
radio emission, and “XDINS” for systems that are isolated NSs with pulsed thermal 
X-ray emission but no detectable radio emission. 


The pulse radiation passes through the ionized interstellar medium (ISM) before 
being detected by radio telescopes. The ISM causes radio signals to disperse in 
the propagation, so the timing of the signal is frequency dependent. For signal of 
frequency f, the time delay t caused by ISM follows t = D x DM/f?, namely that 
the signals of higher frequency arrive earlier than those of lower frequency. The 
parameter D = e?/2am,c = (4148.808 + 0.003) MHz” pce~!cm?s is the dispersion 
constant, where e and me are the charge and mass of an electron respectively. “DM” 
is the dispersion measure, which controls the magnitude of the delay, defined via 
DM = i. nedl, where d is the distance of the pulsar and ne is the electron number 
density. Using the time delay between two different frequencies (fı and f2), 


At =Dx (f° — f} )xDM, (1.1) 


we can calculate the value of DM. With a Galactic electron density distribution 
model, we can calculate the distance of the pulsar. There are two mainstream mod- 
els, the NE2001 model [20] and the YMW model [127]. As post processing for 
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the observation, to rectify the dispersion effect during propagation, one can adopt 
“incoherent dedispersion” or “coherent dedispersion” techniques. 

In reality, the ISM is highly turbulent and inhomogeneous, and the irregularity 
induces extra effects on the propagating signals. It causes multi-path scattering of 
the radio signals, leading pulse profiles of pulsars to have “exponential tails” which 
blur the profile edges and worsen the timing precision. The scattering effect de- 
creases with frequency roughly as f~*, so using a higher frequency band to observe 
reduces the scattering effect [73]. In addition, the relative motion between the pul- 
sar, the scattering screen, and the radio telescope can cause interstellar scintillation. 
It leads to intensity variations of the pulsar’s signals on various timescales and fre- 
quencies [73]. 

After accounting for the dispersion effects of pulse signals, we get a series of sin- 
gle pulses. Figure 1.2 shows 100 continuous single pulses versus spin phase from 
PSR J2222—0137, which was observed by the Five-hundred-meter Aperture Spher- 
ical radio Telescope (FAST). In the figure, the pulse-to-pulse variation of profiles is 
not stable. Generally, the profiles of single pulses of a pulsar are stochastic. But some 
pulsars exhibit unique temporal characteristics, such as sub-pulse drifting [37, 115], 
mode changing [8], nulling [6, | 12], giant pulse [105, 48], and microstructures [21]. 
The investigation of pulsars’ single pulses helps one reveal the radiation mechanism 
of pulsars and the structure of the emission region. Here we do not provide a de- 
tailed introduction on the single-pulse study, more relevant studies on this topic can 
be found in Refs. [95, 70]. We will concentrate on introducing the pulsar timing 
technique that helps us measure pulsars’ physical properties. 

For pulsar timing, we measure the times of arrival (TOAs) of a pulsar’s radio 
signals at radio telescopes on the Earth and use a timing model to fit the TOAs to 
get a phase-connection solution [30, 117]. To generate TOAs, we should use pulse 
profiles to cross correlate with the “standard profile” of the pulsar at that observing 
frequency. To achieve high-precision TOAs, we need to use stable profiles of pulses. 
The shape of profiles is highly variable from pulse to pulse. The integration of hun- 
dreds or thousands of pulse periods leads to stable profiles. In addition, except for 
bright pulsars, single pulses of a pulsar are faint, and it is difficult to detect single 
pulses. So, to get a stable profile with a high signal-to-noise ratio (SNR), we need 
to fold hundreds or thousands of pulse signals according to a pulsar’s spin period. 
In Fig. 1.2, the profile in the top is the integrated profile of PSR J2222—0137, and 
it is folded from the 100 single pulses. Generally, we cannot simply get the spin pe- 
riod directly from the pulse TOA intervals. The propagation effects on signals and 
the motion of pulsars cause various time delays to TOAs. We need to adopt a spe- 
cific timing model to predict these time-delay effects and get the intrinsic rotation 
periods. 

The intrinsic rotation frequency of a pulsar can be Taylor-expanded by [30, 73], 


V(T) = vo + Vo(T —t0) + 5 0(T — to)? +>, (1.2) 


where T is pulsar proper time, Vo = V (to) is the spin frequency at reference epoch 
to, Vo and Vo are the first and second time derivatives of spin frequency at fo. In the 
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Fig. 1.2: Stack of 100 single pulses from PSR J2222—0137, observed with the 
FAST telescope. The line in the top is the integrated profile of these 100 pulses. 


magnetic dipole model, the radiation of a pulsar takes away the rotation energy and 
makes the pulsar spin down, which contributes to Vo and Vo. We can also rewrite 
Eq. (1.2) as the relation between the pulsar’s proper time and pulse number N [30], 


1 ig (T ig)” Eee (1.3) 


Vo (T to)” + 6 


1 
N =No + vo (T to) + 5 


Equations (1.2) and (1.3) describe the pulsar rotation in a reference frame co- 
moving with the pulsar. For us, what we measure is TOAs of pulses recorded by an 
atomic clock on the Earth. The recorded TOAs are in the non-inertial frame of ref- 
erence where the radio telescope is located. So we need to translate the TOAs to an 
inertial frame, and the Solar System Barycentre (SSB) is an excellent approximate 
inertial reference frame. For a solitary pulsar, its proper time T equals to its time at 
SSB of infinite frequency up to a constant. It can be expressed by [30], 


x DM 
f + Ar, +As, + Ag, ; (1.4) 


D 
T = tele tte — 
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where ftele is the time that a pulse arrivals at a telescope, te is the relative correction 
between an atomic clock of the telescope to an average reference clock, D x DM / f? 
is the correction of dispersion effect, f is the observing frequency, Agro is the Römer 
delay caused by the motion of the Earth around the Sun, Ase is the Shapiro delay 
which is a relativistic delay from the curvature of spacetime caused by the presence 
of masses in the Solar system (notably the Sun and the Jupiter), Ag, is the Einstein 
delay which is a combined effect of gravitational redshift caused by masses in the 
Solar system and time dilation caused by the motion of the Earth. 

In gravity tests with radio pulsars, binary pulsar systems are the main objects of 
investigation. For a pulsar in a binary system, the orbital motion of the pulsar and the 
gravity field of the companion lead to extra time delays. For binary pulsar systems, 
Eq. (1.4) is extended to [30], 


T = hele tle — 2x + Ar, +As, +AEo +Arp + Asp +Azgpt+Aap. (1.5) 
Compared with Eq. (1.4), it has four additional terms: Agg is the Römer delay 
caused by the orbital motion of the pulsar; Asp is the Shapiro delay caused by the 
gravity field of the companion star in the binary; Agp is the Einstein delay which 
is caused by the time dilation and gravitational redshift, and they are due to the or- 
bital motion of the pulsar; Aap is the aberration delay that results from aberration 
of the rotating pulsar beam, and this effect depends on the time-varying transverse 
component of the orbital velocity of the pulsar [73]. 

The motion of a pulsar in a binary orbit can be described by six Keplerian pa- 
rameters, and they are the orbital period of the pulsar P,, the orbital eccentricity 
e, the projected semi-major axis x, the longitude of periastron relative to the as- 
cending node @, the time of periastron passage Tọ, and the position angle of the 
ascending node ase which we generally cannot measure solely from pulsar timing 
unless the binary pulsar is very nearby [73]. When the six Keplerian parameters are 
determined, we can describe the motion of a binary system well. However, these 
parameters change with time because of various reasons. For example, for a binary 
pulsar system, the pulsar and the companion are slowly inspiraling and the system 
radiates GWs outward. So we need to consider orbital evolution over time. 

After the first binary pulsar system PSR B1913+16 was detected [56], how to 
provide a simple timing model with correct binary motion became a hot topic that 
was constantly discussed at that time. Considering the first-order post-Newtonian 
(PN) approximation and combining with the pulsar observational data, Damour 
and Deruelle [23, 24] proposed a phenomenological parametrization—the so-called 
parametrized post-Keplerian (PPK) formalism—and later Damour and Taylor [30] 
extended it. The PPK parameters reflect orbital effects beyond the Keplerian or- 
bit that can be extracted from pulsar timing and do not depend on specific boost- 
invariant gravity theories at the 1 PN order. The Damour-Deruelle (DD) timing 
model gets the orbital motion information by fitting a series of Keplerian and PPK 
parameters. The commonly used observable PPK parameters are [24, 30], 


{p5} = {0, y, È, r, s, ôo}, (1.6) 
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where @ is the advance rate of periastron, È, is the decay rate of orbital period, y is 
the amplitude of Agg, r and s are the PPK parameters describing Asg, dg is a PPK 
parameter in App. 

Considering GW radiation makes orbital parameters change with time. Based on 
the Kepler’s equation, Damour and Deruelle [23, 24] provided a Kepler-like equa- 
tion, 


u—esinu = 


(T-T) (1.7) 


2 P, 


P, e 


where u is the eccentric anomaly, and its relation to the true anomaly Ae (u) is, 


1 1 
Ae(u) = 2 arctan y tan J í (1.8) 


The longitude of periastron @ can be modified as, @ = @) + ®Ae(u)/np, where 
np = 27 /P, is the average angular velocity of the orbit. 

Damour and Deruelle [24] provided the time delay terms for binary systems. The 
Römer delay is, 


1/2 
Arg = xsin@ [cosu — e (1 + ô,)] +x [1 -e° (1 +80)? cos@sinu, (1.9) 


where ô, is also a PPK parameter, but we cannot get it by directly fitting it in the 
timing model, because it can be absorbed by redefining other parameters. dg and 
6, describe the relativistic deformations of the orbit. Damour and Deruelle used 
three eccentricities, and the two new eccentricities are defined via e, = e (1 + ô) 
and eg =e (1 + ôo). 

The Einstein delay of a binary is, 


Agg = Ysinu, (1.10) 


where y represents the amplitude combining time dilation and gravitational redshift. 
The Einstein delay is degenerate with the Romer delay, and the degree of degeneracy 
depends on the change of @ [30], namely, if @ changes significantly, the degeneracy 
can be broken. 

The Shapiro delay of a binary is [12], 


Asp = —2rin{] —ecosu—s [sin w(cosu—e) + (1 - e)" ? cos osinu| } , (1.11) 


where r and s are the above-mentioned range and shape parameters respectively. 
The Shapiro delay effect will be noticeable when the orbital inclination i, which is 
the inclination of the orbit with respect to the line of sight, is close to 90°, that is, the 
binary system is close to edge-on. Otherwise, this effect can be partially absorbed 
in other timing parameters. A strong Shapiro delay effect provides a direct measure- 
ment of the mass of the companion and the inclination angle 7. For low-eccentricity 
binary systems, the PPK parameters r and s have a high correlation, especially when 
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the inclination angle / is not sufficiently close to 90°. In this situation, we cannot get 
rand s with high precision. Freire and Wex [41] provided two new PPK parameters, 
hz and h4, to alleviate this problem. The PPK parameters h3 and hy are the ampli- 
tudes of the third and fourth harmonics from the Shapiro delay’s Fourier expansion. 
Compared with r and s, h3 and h4 are less correlated with each other. For edge-on 
systems, parameters h3 and € are introduced to replace h3 and h4, where € = h3 /h4 
is the ratio of amplitudes of successive harmonics. The parameters h3 and ¢ provide 
a superior description of the constraints on the orbital inclination i and the masses 
of the binary than r and s. In the new parametrization, the Shapiro delay is rewritten 
as [41], 

Asp =—Fyln(1+6?7—2¢sin®) , (1.12) 
where @ is the longitude relative to the ascending node. The relations between 
(h3, Ç) and (s, r) are s =2€/(€7 +1) andr = h3 /¢? [41]. 

Finally, the aberration delay is [73], 


Aap =A {sin [@ +A. (u)] +esin@}+B{cos[@+A,(u)|+ecos@}, (1.13) 


where A and B are both PPK parameters, but they are not separately measurable by 
pulsar timing. 

After considering these time delay terms contributed to the timing model, we can 
precisely measure orbital parameters by pulsar timing. For boost-invariant gravity 
theories, in the DD timing model, the PPK parameters are functions of Keplerian 
parameters, masses of the binary and parameters of gravity theories. In GR, PPK 
parameters are functions of Keplerian parameters and masses of the binary only. In 
GR, PPK parameters in Eq. (1.6) are expressed as [73, 117], 


; 3np v? 
o= 5, (1.14) 
N me VŽ 
y=< (1 me) me 7 (1.15) 
np m/ mc 
a =. (1.16) 
E 
s=sini =n E, (1.17) 
Me Vp 
7/2)m2 +6m,me +2m2 v2 
j- a ae (1.18) 
m E 
, 1927 m pm... Vp 
ae caer ee (1.19) 
m C 


where V, = (Gmnp)'/? and f(e) = [1 + (73/24)e? + (37/96)e*] (1 — e)~7/?. In the 
above equations, mp is the mass of the pulsar, mc is the mass of the companion, m 
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is the total mass of the binary and x = ap sini/c is the projected semi-major axis of 
pulsar’s orbit in the unit of seconds. 

In relativistic gravity theories, the spins of the bodies that orbit in a binary sys- 
tem can affect their orbital and spin dynamics. Three contributions are induced by a 
rotating body moving in a binary pulsar system at the leading orders. They are the 
spin-orbit interaction between the spin of pulsar S, and the orbital angular momen- 
tum L, the spin-orbit interaction between the spin of companion S, and L, and the 
spin-spin interaction between S, and S, [117]. The spin-spin interaction is generally 
many orders of magnitude smaller than the 2 PN terms and the spin-orbit effects, so 
that one can ignore it. We only introduce the spin-orbit interaction here. In a bi- 
nary pulsar system, the spin-orbit interaction leads the orbit and the two spins to 
precess [117, 7]. The main effect of the precession of the spin is a secular change 
in the orientation of the spin, which is caused by the effect of spacetime curva- 
ture and the precessing rate is independent of the spin. This effect is the so-called 
geodetic precession, and it has been detected in some binary pulsar systems, such 
as PSR JO737—3039 [15] and PSR J1906+0746 [32]. 

For a pulsar in a binary system, geodetic precession causes the pulsar’s spin to 
precess around the total angular momentum. As the orbital angular momentum is 
normally much larger than the pulsar’s spin, we can regard that the pulsar’s spin 
precesses around the orbital angular momentum. In GR, the average rate of geodetic 
precession of the pulsar’s spin at the leading order is [73], 


c VZ 
TEK (e) b (1.20) 


~ ]=e 2mp m? e 


The secular change of the spin orientation causes changes in the observed emission 
properties of the pulsar because the line-of-sight will cut through different regions 
of the magnetosphere. One can detect a changing integrated profile and polarization 
angle (PA) in data. In turn, one can use the change of profile or PA to calculate Qp, 
and compare the results with GR’s predictions [32]. 

The orbital precession has two effects that need to be considered. The first effect 
is due to the mass-mass interaction, which belongs to the PN terms of two point 
masses. The second effect is due to the Lense-Thirring effect related to the spins of 
two stars [7]. Damour and Schaefer [29] decomposed the precession of the orbit in 
terms of observable parameters of pulsar timing. One of the parameters is @,, which 
can be decomposed into, 


Op = OPN + ÖLT, + ÒLT, , (1.21) 


where Wpn is the PN term and the 1 PN order is given in Eq. (1.14); OT, and OLT, 
are due to the Lense-Thirring effect of the pulsar and the companion respectively. 
For a binary NS system, the spin of the companion is generally slow, so the @rr, 
that comes from the contribution of the companion’s spin can be ignored. If S, is 
parallel to the orbital angular momentum L, the @pr, can be simplified to [7, 54], 
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3 
. = 3ny SoS V; 
OLT, = -7 2PpS8pll 23? 


(1.22) 


—1/2 : 
5 / , and J, is the 


moment of inertia of the pulsar which depends on the equation of state (EOS) of 
NSs. As we can see, the value of OLT, depends on J,. So if we can measure OLT,» 
we can get a value of J, which can help us determine the EOS of NSs. However, 
using the best-measured double pulsar as an example [61, 54, 60], the value of 
OT, is at the same order as the value of @pn, so only highly relativistic binary 
pulsar systems with high-precision timing results might have the ability to limit the 
moment of inertia of NSs with this method. 


where X, = mp/m, Bs = clpnp/Gm,, gp = (3X? +Xp) (1 -e) 


1.3 Hulse-Taylor Pulsar, Double Pulsar, and Triple Pulsar 


With the theoretical background in the last section, we now introduce three of the 
most famous radio pulsars in the field, the Hulse-Taylor pulsar PSR B1913+16 [56, 
109, 114], the double pulsar PSR JO737—3039A/B [16, 76, 59, 15, 60], and the 
triple pulsar PSR JO337+1715 [87, 5, 111]. 


1.3.1 The first binary pulsar: PSR B1913+16 


PSR B1913+16, a recycled pulsar in a highly eccentric 7.75-hr orbit, was the first 
pulsar detected in a binary system. The binary is a double NS system, and it was 
discovered using the Arecibo telescope in 1974 [56]. We also call this pulsar the 
Hulse-Taylor pulsar in honor of its discoverers. The discovery of this system in- 
directly proved the existence of GW radiation for the first time [109]. The latest 
timing results of this pulsar are published in Ref. [114], where 9257 TOAs were 
acquired over 35 years, and the analysis provided stringent tests of gravity theo- 
ries. The orbital parameters in DD timing model are listed in Table 1.1. Three PPK 
parameters—the advance rate of periastron @, the amplitude of Einstein delay y, and 
the decay rate of orbital period P,—were well determined [114]. In addition, two 
PPK parameters, r and s, of the Shapiro delay effect from this system were mea- 
sured. All these measurements are consistent with GR predictions. The relativistic 
shape correction to the elliptical orbit, oe is also detected for the first time. 
Assuming GR is the true theory of gravity, Eqs. (1.14—1.15) were used to get the 
masses of the binary, mp = 1.438 + 0.001 Mọ and m; = 1.390 + 0.001 Mo [114]. 
Inserting the measured Keplerian parameters and the derived masses into Eq. (1.19), 
Weisberg and Huang [114] got a theoretical prediction in GR, por = — (2.40263 + 
0.00005) x 107!?. Table 1.1 shows the value of the observed decay rate of orbital 
period pobs, and it must be corrected by terms ponk, the dynamical contributions 
from the Shklovskii effect and poal the differential Galactic acceleration between 


12 Hu et al. 


Table 1.1: Orbital parameters of PSR B1913+16 [1 14]. 


Parameter Value 

x= ap sini/c (s) 2.341776(2) 

e 0.6171340(4) 

P, (d) 0.322997448918(3) 
@ (deg) 292.54450(8) 

@ (degyr—!) 4.226585(4) 

y (ms) 0.004307(4) 

pos —2.423(1) x 107!? 
obs 4.0(25) x 1076 

s 0.6870 10 

r (us) 9.6735 

4 0.38(4) 

hz (s) 0.6(1) x 1076 


the SSB and the pulsar system. After subtracting the total contribution from pe and 
PS — (0.025 + 0.004) x 10-!, an intrinsic value is obtained, Pi" = —(2.398 + 
0.004) x 107! [114]. Using pR and pintr, Weisberg and Huang found that the 
ratio of the observed orbital period decrease caused by the GW damping to its GR 
prediction is [114], 


pintr — (2.398 + 0.004) x 107!? 
B = 0.9983 + 0.0016. 1.23 
PSR  — (2.40263 + 0.00005) x 107!2 oe 


This result provided the most precise test of GW emission in 2016. One can no- 
tice that the uncertainty of Pi" is dominated by the errors of Pe** and PS, which 
depend on the precision of the pulsar distance measurement and the Galactic accel- 
eration model. 

As mentioned above, assuming GR is correct and using the measured two PPK 
parameters, one can get mp and me. If one measures three or more PPK parame- 
ters, one can do self-consistency checks in GR. Mass-mass diagram illustrates the 
self-consistency checks. Using Eqs. (1.14-1.19), one can plot the m,-m, curves for 
each observed PPK parameter and check whether they intersect at a common re- 
gion within measurement uncertainties. Figure 1.3 shows the mass-mass diagram of 
PSR B1913+16 [114], and the curves of five PPK parameters all meet within mea- 
surement uncertainties. For PSR B1913+16, three PPK parameters, @, y and ae 
are measured with high precision; the 1-o errors of them are smaller than the widths 
of the curves in the figure. All lines intersect in a tiny region, which provides a strict 
test of GR in the strong-field condition. 
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Fig. 1.3: Mass-mass diagram of PSR B1913+16 [114], assuming that GR is right and 
using the five PPK parameters listed in Table 1.1. The width of each band represents 
the +1-o uncertainty. For ©, y and por. the uncertainties of them are smaller than 
the widths of the lines. The intersection of all bands is consistent within a small 
region, which illustrates the agreement of the observations with GR. 


1.3.2 The double pulsar system: PSR J0737—3039A/B 


PSR J0737—3039A/B, discovered in 2003, is currently the only known double pul- 
sar system in which both NSs are detectable as radio pulsars [16]. PSR JO737—3039A 
is a 22.7-ms MSP pulsar which was first-born and spun up by accreting materials 
from the companion; PSR JO737—3039B is the second-born pulsar with P = 2.7 s 
and, due to the geodetic precession, its radio emission disappeared in 2008 [15]. The 
two pulsars orbit each other in a mildly eccentric (e = 0.088) orbit with P, = 2.45 hr 
[76, 59]. The double pulsar system’s orbital period is only one-third of that of the 
Hulse-Taylor pulsar, which means PSR JO737—3039A/B has a higher average or- 
bital velocity and a larger acceleration. So it is a more relativistic binary system and 
an excellent testbed for strong-field gravity when well timed. 

In 2006, Kramer et al. [59] reported 2.5-yr timing results of PSR JO737—3039A/B, 
which provided four independent strong-field tests of GR. Although only using 2.5 
years of timing data, because of the larger relativistic effect in this system, the PPK 
parameters ©, y and P, were already well measured. The orbit is nearly edge-on, 
so the PPK parameters r and s of the Shapiro delay effect also have precise mea- 
surements. Because the two pulsars’ emissions are detectable for this system, the 
projected semi-major orbital axes xq and xg are both measured with high preci- 
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Table 1.2: Four independent tests of GR provided by the PPK parameters of the 
double pulsar [59]. The last column shows the ratios of observed values to their 
expectation in GR. 


PPK Parameter Observed Expected in GR Ratio 

È, (107!) —1.252(17) —1.24787(13) 1.003(14) 

y (ms) 0.3856(26) 0.38418(22) 1.0036(68) 
s 0.99974(—39, +16) 0.99987(—48, +13) 0.99987(50) 
r (us) 6.21(33) 6.153(26) 1.009(55) 


sion. Therefore this system provides an accurate measurement of the mass ratio 
R = ma/mp = xg /xa, which is a theory-independent parameter in boost-invariant 
gravity theories. Kramer et al. [59] used the most precisely measured PPK parameter 
@ = 16.89947(68)degyr—! and the theory-independent parameter R = 1.0714(11) 
to derive the masses of the binary, ma = 1.3381(7) Mj and mg = 1.2489(7) Mo. 
The masses are then used to calculate the theoretical values of other PPK param- 
eters with the assumption that GR is right. In Table 1.2 we show the results of 
four independent gravity tests. The fourth column is the ratio of observed values to 
their expectation in GR, wherein the ratio for s was the best available Shapiro-delay 
test of GR in the strong-field region at that time. With only 2.5-yr data, it reveals 
PSR J0737—3039A/B’s outstanding ability to test GR. 

For PSR J0737—3039A/B, it shows a roughly 30-s eclipse when pulsar A passes 
behind pulsar B, namely, pulsar A is at the superior conjunction of its orbit. Pul- 
sar A’s pulse signals will be absorbed by the magnetosphere of pulsar B when the 
eclipse happens. For pulsar B, the misalignment angle between the spin vector and 
the total angular momentum vector is significant, dg ~ 50°, leading to a measurable 
geodetic precession of pulsar B’s spin, Sg. The geodetic precession of Sg causes a 
changing absorption effect at the eclipse. Breton et al. [15] used a simple geometric 
model to characterize the observed changing eclipse morphology and measured the 
geodetic precession of Sg around the total orbital angular momentum. They got a 
relativistic spin precession rate, Qg = 4.77°(+0.66°, —0.65°) yr}, and this result 
is consistent with the value, aoe = 5.0734(7)° yr—!, predicted by GR, providing a 
new test of GR in the strong-field regime [15]. 

In 2021, Kramer et al. [60] reported 16.2-yr timing data of the double pulsar 
observed with six radio telescopes. These results show an excellent ability to test 
gravity theories because they not only provide very good test of GR in the strong 
field but also reveal new relativistic effects that are observed for the first time. These 
timing results provide seven PPK parameters, which are more than those from any 
other known binary pulsars, and some of the parameters need to take higher-order 
contributions into account. Our following discussion will highlight some novel as- 
pects of the double pulsar from this new set of data. 
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As mentioned before, the DD model is a compact timing model that uses the 
solution of the 1 PN two-body problem. However, the precision of timing for the 
double pulsar is so high that higher-order corrections to the timing model need to be 
added. For @, 2 PN and the Lense-Thirring effects are needed [54] 


O = Oi pN + @ pnt ÖLT, ;» (1.24) 


where @pn = @1pnfoV?/c? with 


1 (3 3. 27 5 23 1 
fo=i— Ge paat n) (gx XA i) (1.25) 


where Xa = ma/m. For PSR JO737—3039A, the misalignment angle of S4 relative 
to the total angular momentum vector is very small, ôa < 3.2° [38]. So Eq. (1.22) is 
used to calculate @pr,. Using values from timing, it was found that @2pn ~ 4.39 x 
10-4 degyr—! in GR, about 35 times of the measurement error of @ (see Table 1.3). 
For @r,, except that the moment of inertia of pulsar A (74) depends on the EOS of 
NSs, all quantities in Eq. (1.22) are measured with high precision, and we get [60], 


Orr, ~ 3.77 x 1074 x 1) deg yr!, (1.26) 
where 1) = I,/(10* gcm?). Generally, 1) is of order unit, so ®t, and @2pn 
have the same order of magnitude [61]. 

In Eq. (1.9), the Römer delay depends on the projected semi-major axis x. In 
Newtonian gravity, the binary mass function describes the relation between the in- 
clination angle i and x by Eq. (1.17). Because of the high precision measurements, 
the mass function now needs to extend to 1 PN [60], 


NpxX C 


sini = 14/3 1x5, \% (1.27) 
aaa (a la qe ye |? 


with Xg = mp/m. Such an extension was the first time for any binary pulsar system. 

GW damping enters the GR equations of motion at 2.5 PN level [c.f. Eq. (1.19)]. 
For the double pulsar, it is extended to 3.5PN, with P3°PN = p25PNX3.5PN ~ 
—1.75 x 10~!7; for the factor X3°P% one can refer to Eqs. (12-13) in Ref. [54]. 
This value is about 4.5 times smaller than the timing precision of P;,. The mass loss 
caused by the spin-down of the pulsars can also contribute to the observed P,, and 
for pulsar A, it is te =2.3x 1071 x 1), Using constraints from Ref. [33] and 
the radius-J, relation in Ref. [65], one gets pra = 2.9(2) x 10717, which is 3 times 
smaller than the error of Pi". Besides Pi", the non-intrinsic contributions PS" and 
poal also need to be considered. These two terms are both functions of the pulsar 
distance d. The parallax from VLBI and pulsar timing was combined to provide a 
weighted distance, d = (735 + 60) pc. Considering all the above, P, is expressed as, 


P, = pees pen + pa + Pe, (1.28) 
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Table 1.3: Orbital parameters of PSR JO737—3039A [60]. 


Parameter Value 

Projected semi-major axis, x (s) 1.415028603(92) 
Eccentricity, e 0.087777023(61) 
Orbital period, P, (d) 0.10225 15592973(10) 
Epoch of periastron, Tọ (MJD) 55700.233017540(13) 
Longitude of periastron, @p (deg) 204.753686(47) 
Periastron advance rate, @ (deg yr!) 16.899323(13) 
Einstein delay amplitude, y(ms) 0.384045(94) 

Change rate of orbital period, P, —1.247920(78) x 107!? 
Logarithmic Shapiro shape, Zs 9.65(15) 

Range of Shapiro delay, r (us) 6.162(21) 

NLO factor for signal propagation, gnLo 1.15(13) 

Relativistic deformation of orbit, dg 13(13) x 1076 
Change rate of projected semi-major axis, x 8(7) x 10716 

Change rate of eccentricity, ¢(s~!) 3(6) x 10716 


Derived parameters 


sini = 1 — exp(—zs) 0.999936(+9, —10) 
89.35(5) or 90.65(5) 
2.587052(+9, —7) 
1.338185(+12, —14) 


1.248868(+13, —11) 


Orbital inclination, i (deg) 
Total mass, m (Mo) 
Mass of pulsar A, ma (Mo) 


Mass of pulsar B, mp (Mo) 


where Pe*t = PSK 4 pGal — _1 68 ttl x 10716, leaving the intrinsic contribution, 
pint = P, — Pext = —1,247752(79) x 107!? [60]. The error of Fi" is still dominated 
by the error of P,. In the future improving the timing precision of the observed P, 
will improve this gravity test. 

The Shapiro delay effect in Eq. (1.11) is the leading-order effect caused by the 
companion’s mass influencing the signal propagation. It was obtained by integrat- 
ing along a straight line and assuming a static mass distribution [12]. For the double 
pulsar one needs a lensing correction to the Shapiro delay. It restores the fact that 
pulsar A’s radio signal propagates along a curved path. Equation (1.11) is rewritten 
as Asg = —2rInA,, and Kramer et al. [60] added a term corresponding to the lens- 
ing correction 8A] = 2rc/apr, where ag is the semi-major axis of the relative orbit. 
Considering a non-static mass distribution, a 1.5 PN correction oA from the re- 
tardation effect was taken into account. So the signal propagation delay is extended 
to Asg = —2rln (Ay +84 + 5AM). For the aberration delay Aap, considering 
the rotational deflection delay, Kramer et al. [60] used a lensing correction term, 
Dcos P/A, where ® = œ + A,(u). Nevertheless, the next-to-leading-order (NLO) 
contributions from the Shapiro delay and the aberration delay cannot be tested sepa- 
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rately in the double pulsar system. These contributions are rescaled with a common 
factor quro to test the significance of NLO signal propagation contributions. The 
factor quro can be fitted in the extended timing model with the double pulsar data. 

After providing an extended timing model with higher-order contributions, one 
can use the timing results in Table 1.3 to determine the masses and test GR. Gener- 
ally, one needs two PPK parameters to determine the masses of a binary system in 
GR. However, when including higher-order contributions, the parameters @ and a 
are also functions of I4, namely @ = @(m,,mp,1,) and Pi" = P!™ (m4 ,mp,I,). Us- 
ing constraints on the EOS obtained from the binary NS merger event GW170817 
[1], one has restriction on 74. Kramer et al. [60] used the information of 74 and the 
timing values of @ and s and got the masses of the binary, 


ma =1.338185(+12, —14)Mo, (1.29) 
mp = 1.248868(+13, —11)Mo, (1.30) 
m =2.587052(+9, —7) Mo . (1.31) 


If one ignores any existing constraints on the EOS of NSs and simultaneously 
determine m4, mg and I4, one needs three PPK parameters. With ©(m4,mp, I4), 
Pi"t(m4,mpg, Ia) and s(m,,mg), one not only gets the masses of the system, but also 
constrains Z4 to Z4 < 3.0 x 10% gcm? with 90% confidence, complementing the con- 
straints from LIGO/Virgo [1, 64] and NICER [103] observations. 

The double pulsar provides seven PPK parameters, {@, y, Py, r, s, Qg, 59}, and 
a theory-independent parameter R (c.f. Table 1.3). For dg, its measured value is not 
significant, so it is not used in testing GR. When a value for Z4 is chosen, these 
PPK parameters are functions of binary masses and those well-measured Keplerian 
parameters. So using six measured PPK parameters and the mass ratio R, one can 
perform five independent tests of GR, as shown in Fig. 1.4 [58]. 

As mentioned earlier, including the 3.5 PN contribution from GW damping and 
removing the mass loss and non-intrinsic contributions, one has [60], 


PGW — pint _ pra — —1,247782(79) x 10-™. (1.32) 
Using the masses shown in Eqs. (1.29-1.30), the predicted value from GR is, 


pGW.GR _ pGW.GR(2.5PN) 


È 5GW,GR(3.5PN) _ 
b b 


+F, = —1.247827(+6, —7) x 10°. 
This value provides the most precise test for the GW quadrupolar emission [60], 
Pew JÈS WOR _ 0,999963 (63). It is also the most precise result among the different 
independent tests that one can obtain from this system. It is in agreement with GR 
at a level of 1.3 x 1074 with 95% confidence [60]. Other independent tests are listed 
in Table 1.4. PSR J0737—3039A/B provides the strictest limit on GR than any other 
binary pulsar systems. 
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Fig. 1.4: Mass-mass diagram of PSR JO737—3039 [60]. It shows the constraints 
from six PPK parameters and the mass ratio R. The inset shows the zoom-in region 
of the intersection of these measurements. The solid black line for @ responds to 
i) = 1.32. The gray band indicates the range for @ with a changing ee The 
dashed black line for @p does not contain the contribution of the Lense-Thirring 


effect. 


1.3.3 The triple system: PSR J0337+1715 


Binary pulsar systems provide us with great opportunities to test GR. There are also 
some pulsars that exist in more complicated systems. For example, PSR B1257+12, 
a pulsar moves in a system that contains at least three low-mass planets. This pul- 
sar helped people discover the first exoplanet system [121, 77]. Here we focus on 
PSR J0337+1715, which is a 2.7-ms pulsar in a hierarchical triple system with two 
white dwarf (WD) companions [87]. The inner pulsar-WD system is a 1.6-d circular 
orbit system and the mass of the WD is 0.19751(15) Mo. The outer orbit is wider, 
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Table 1.4: Relativistic effects measured in the double pulsar system and the resulting 
independent strong-field tests of GR [60]. The right column shows the ratio between 
the observed quantity and its prediction in GR. 


Parameter Ratio 
Shapiro delay shape, s 1.00009(18) 
Shapiro delay range, r 1.0016(34) 
Time dilation, y 1.00012(25) 
Periastron advance, ®© 1.000015(26) 
Change rate of orbital period, Ê, 0.999963(63) 
Orbital deformation, ôg 1.3(13) 

Spin precession, Qg 0.94(13) 


Tests of higher-order contributions 


Lense-Thirring contribution to @, Aur 0.7(9) 
NLO signal propagation, qno [total] 1.15(13) 


with P, = 327.2d, and the WD’s mass is 0.4101(3) Mo [87]. The special system 
composition and high-precision timing results of PSR JO337+1715 make it an ideal 
laboratory that can provide a strong equivalence principle (EP) test of GR [5, 111]. 

The principle of equivalence has played an important role in the development of 
theories of gravity. In Newton’s gravity, the principle of equivalence is the corner- 
stone of the theory. It assumes that the inertial mass equals the passive gravitational 
mass. Now, this is so-called the weak EP (WEP), and its alternative statement is 
that the trajectory of a freely falling test body is independent of the body’s internal 
structure and composition. With the development of gravity theory, the principle of 
equivalence has been extended. Except for satisfying the WEP, one believes that a 
gravity theory should also satisfy the local Lorentz invariance (LLI) and local posi- 
tion invariance (LPI). Combined with WEP, we call it the Einstein EP (EEP). In the 
EEP, LLI states that a non-gravitational experiment does not depend on the velocity 
of the freely falling reference frame, and LPI states that the outcome of any local 
non-gravitational experiment is independent of when and where the experiment is 
performed. Metric theories of gravity satisfy the EEP [118]. The strictest bound 
comes from the strong EP (SEP), it demands that WEP is valid for self-gravitating 
bodies as well, and a gravitational experiment needs to embody LLI and LPI simul- 
taneously. For now, GR seems to be the only viable metric theory that satisfies SEP 
completely. So testing SEP is equal to testing GR in this sense. 

There are many methods to test SEP, and one of them is considering the uni- 
versality of free fall of two self-gravitating bodies in an external gravitational field. 
Some tests from the Solar System have been done and provided tight constraints 
on the violation of the SEP, but the bodies in the Solar System are weak in terms 
of self-gravity. For instance, the fractional binding energy difference between the 
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Earth and the Moon is ~ —4 x 107!°. Pulsar-WD systems can provide an evident 
fractional binding energy difference between the NS and the WD, at the order of 
—0.1, which is many orders of magnitude larger than the Earth-Moon system’s dif- 
ference. If the acceleration difference between the pulsar and the WD is present, 
it will lead to a characteristic change in its orbital eccentricity evolution. However, 
a pulsar-WD system’s free fall in the Galactic gravitational field cannot provide a 
tight limit on the SEP violation parameter A, because the Galactic gravitational ac- 
celeration at the location of the pulsar is comparably weak, ~ 2 x 1078 cms~?. For 
example, PSR J1713+0747 provides a constraint, A < 0.002 [133]. For the triple 
system PSR JO337+1715, the outer WD can provide an external gravitational accel- 
eration for the inner PSR-WD system, and its strength is ~ 0.17cms~2, which is 
~ 10’ times larger than the Galactic gravitational acceleration. So PSR J0337+1715 
makes a significant improvement on the constraint of the SEP violation parameter, 
A < 2.6 x 1076 [5] and A = (0.5+0.9) x 10~ [111] by two independent groups. 


1.4 Neutron Star Structures 


As introduced in the previous section, pulsars, which are rotating NSs, can provide 
unique possibilities to probe the strong-field region of GR. With a typical mass of 
1.4 Mo and a radius of 12 km, NSs not only provide the strong-field environment of 
gravity, but also are the most extreme laboratory of particle physics. In this section, 
we will introduce the basic theory describing the NS structures, and show some 
properties of a single NS that can be measured with pulsar timing technique and 
other observations. Differing from the tests of gravity from binary pulsar systems, 
which use the dynamic aspect of the gravity theory, testing gravity with a single NS 
provides the properties of an equilibrium state. In the following, we take G = c = 1 
except when the units are written out explicitly, and the convention of the metric is 
(—,+,+,+). 

The basic information of a compact star is its mass and radius, which reveal the 
inner structure and interaction among the star’s components. In GR, to derive the 
structure of a spherical, non-rotating NS, we begin with the Einstein’s field equa- 
tions, Guy = 827 ,y, together with some specific matter model, for example a perfect 
fluid, TY” = (p +p)u"u’ + pg!” . Then one adopts the metric of a static, spherically 
symmetric spacetime 


2 dr? 
i. d? — T prado, (1.33) 
r 1—2m/r 


and inserts it into the field equations. The equation of a star is described by the 
Tolman-Oppenheimer-Volkov (TOV) equation [84, 110] 


dp (p+p)(m+4ar*p) 
dr r(r—2m) , ae) 
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where m(r) is called the mass inside the sphere of radius r, defined by the inte- 
gral, m(r) = fj 47r’pdr, just as in the Newtonian theory. Results show that for any 
physical realizable matter, m(r)/r < 0.485, so that the right hand side of Eq. (1.34) 
is never singular [13]. The dimensionless ratio between a star’s mass and radius, 
C = M/R, is called the compactness of the star, which represents how compact a 
star is. For Schwarzschild black holes (BHs) in GR this number is equal to 0.5 and 
for a typical NS, the compactness can reach 0.2. As a comparison, the compactness 
of the Sun is only 2 x 1076. This is crucial for the SEP in Sec. 1.3.3. 

To complete the above equations, a relation called the EOS between the pressure 
p and the energy density p is needed. For a simple fluid in local thermodynamic 
equilibrium, there always exists a relation of the form, p = p(p,S), where S is the 
specific entropy. For old NSs that are sufficiently cold or under the assumption of 
constant entropy, we can restrict this relation to a one-parameter functional, p = 
P(p). 

For a given EOS, one can then solve the TOV equation with initial condition 
P|r=0 = Po to get the structure of a NS. With the uniqueness theorems for ordinary 
differential equations, the solutions will form a one-parameter sequence with the 
parameter being the central density. Correspondingly, for a given theory of gravity, 
this sequence, the so-called mass-radius relation, also uniquely determines the un- 
derlying EOS [68]. The total mass of the NS is m(R), where R is the radius of the 
surface defined by p(R) = 0. For a normal NS that is gravitationally self-bounded, 
the density also approaches zero at the star’s surface, i.e. p (R) = 0. But if a NS 
consists of quarks or strangeons that are not gravitationally self-bounded, there can 
be a density discontinuity at the surface of the star [28, 89, 43]. 

In Fig. 1.5 we show some EOSs and the corresponding M-R relations for NSs. 
The EOSs inside the NSs are still not clear today. The main components of NSs are 
mainly regarded as neutrons, but there also should be protons, electrons, muons, and 
so on. It remains an unsolved problem in quantum chromodynamics for calculating 
the super-dense nuclear matter at low temperatures. Furthermore, there are some 
arguments that the ingredients of NSs should be quarks or strangeons [120, 122, 63, 
43]. From Fig. 1.5, we can see that there is large difference between different EOS 
models, especially between the models for a normal NS and models for a quark or a 
strangeon star. For those NSs that are not gravitationally self-bounded, the behavior 
at the low mass limit is like a star with constant density. This feature still holds 
approximately when the mass of the star increases, and leads to various differences 
in other aspects like the moment of inertia and the star oscillation frequency (see 
e.g. Ref. [67]). But one can notice that all of these EOS models have a maximum 
mass in their M-R relations. For WDs this maximum mass is famously known as the 
Chandrasekhar limit, which is around 1.4Mọo. Due to the uncertainty in the EOS, the 
maximum mass of NSs can vary largely, from about 1.5Mo to over 3M,. Known 
as the turning point theory [47], a NS with mass and radius crossed the maximum 
mass point in the M-R relation becomes unstable. So the discover of supermassive 
NSs can help us constrain the EOS of NSs. The two horizontal lines in Fig. 1.5b are 
the observed masses of PSR J0348+0432 [4] and PSR J1614—2203 [39], and then 
EOSs which cannot support a maximum mass higher than 2M, are excluded by 
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Fig. 1.5: EOSs for NSs [66] and the related M-R relations (figure courtesy of Y. 
Gao). The horizontal lines in Fig. 1.5b are the observed masses for two massive 
pulsars [4, 39]. The SQM model is related to quark stars [66]. Curves passing the 
highest point in their M-R relations are cut off, as they represent unstable NSs. 


these observations. But if one considers gravity theories different from GR, the M-R 
relation of the same EOS can vary largely from their GR counterpart. Thus a precise 
measurement of the M-R relation can constrain not only the EOS but also the gravity 
theory [102]. Previous measurements from the NICER satellite by modeling the 
thermal x-ray waveform of the isolate pulsar PSR J0030+0451 give a constraint on 
the mass and equatorial radius of the star simultaneously [90, 81]. But measurements 
from more NSs with higher precision are still needed. 
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The above discussion is based on the assumption of a static star. For real pul- 
sars, they can have spin periods from several milliseconds to tens of seconds (see 
Fig. 1.1). Some studies have shown that fast-rotating NSs can support a larger max- 
imum mass and it will have important influence on the evolution of the binary NS 
merger remnant [9, 132]. One can estimate the maximum angular velocity of a NS 
from the Newtonian theory, Qnax = (M/ R?)!/ 2 which is related to the mass shed- 
ding limit. For a typical NS, this corresponds to a spin period of about 0.5 ms. Thus 
for most pulsars, the approximation of slow rotation is sufficient, and in this case, 
the mass-radius relation of NSs will not change, to the first order in Q. 

The basic property related to rotation is the moment of inertia. In GR, the mo- 
ment of inertia of a slowly rotating NS can be calculated through the slow rotation 
approximation [49]. Under this approximation the metric is writen as 


2 dr? 
ds? = guydx"dx” = — (1 = =~) di? + —__ +r (d0? + sin” 46?) 
r m/r 
+r sin? 0 (dọ + (@(r)—Q)dt) , (1.35) 


where m/(r) is just the same as that solved from the TOV equation for a non-rotating 
NS, Q is the angular velocity of the star and @(r) shows the effect of frame drag- 
ging. An observer falling freely from infinity to r will have an angular velocity 
Q — @(r) to the first order in Q, thus seeing the fluid element at r has an angular 
velocity @(r). Function @(r) is determined by the linear differential equation, 


d? 4 d 4d 
o mr (A) ) o (1.36) 


= | 4w : 
dr2 roan p) ¢ dr zs r dr 


The boundary condition is @ = Q when r goes to infinity. Besides that, in order 
to keep d?/dr* to be finite at the origin, we have d@/dr|,-o = 0. The angular 
momentum J carried by the star can be read out from the asymptotic behavior of gro 
at r+ œ, go = (© — Q )r? sin? O + —2J sin? 0 /r, which indicates that 


_ ido, 


= é 1.37 
6 dr roo ( ) 
The moment of inertia of the star thus is defined to be 
J do 
n a a * 1.38 
Q 6@ dr |r ( ) 


Combining with Eq. (1.36), the moment of inertia of the star can be calculated from 


di(r) _ 24ar*(p +p) 51I P 
dr 3 1—2m/r 2r ro 


(1.39) 


with initial condition /(r = 0) = 0. This equation, similar to the TOV equation, 
can be regarded as the corresponding equation in the Newtonian theory with some 
correction terms that come from GR. 
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Fig. 1.6: Moment of inertia of NSs with different EOSs (figure courtesy of Y. Gao). 
Unstable solutions are cut off as in Fig. 1.5b. 


In Fig. 1.6 we show the NS moment of inertia with some EOSs. It is relatively 
hard to measure the NS moment of inertia directly, because there is no leading or- 
der effect caused by the current in Newtonian gravity. Some indirect measurements 
come from modeling the X-ray profile of X-ray pulsars and combining with proper 
universal relations [103]. Notably, as explained in Sec. 1.3.2, the observations of the 
double pulsar binary starts to constrain the moment of inertia through the spin-orbit 
coupling [60]. Measuring the moment of inertia independently is valuable in gravity 
tests as it relates to the universal relationship that will be discussed later. 

There are some other quantities that describe the NS structure and can be mea- 
sured today. One example that will be used later is the tidal Love number of 
NSs [51]. This quantity describes how easily an object can be tidally deformed and 
it can be measured from the GW signal of colliding NSs [1]. At the early stage of 
the inspiral, the motion of the binary system is dominated by the point-mass dy- 
namics, which is similar to, but with higher-order corrections, what we discussed in 
Sec. 1.2. Near the end of the inspiral, the non-uniform gravity field of the compan- 
ion will cause both stars to be tidally deformed, and the internal degrees of freedom 
of the stars begin to influence the GW signal. The observation of the famous event 
GW170817 set an upper limit on the NS tidal Love number that A < 800, which 
excludes some EOSs that are too stiff [1]. 

There is a problem that comes up when we want to test strong-field gravity with 
NS structures: the uncertainties in the EOS [96, 102]. Nuclear physics in the labo- 
ratory cannot tell us which EOS is correct yet. Thus, for example, the mass-radius 
relation of NSs depends not only on the theories of gravity but also on the choice of 
EOSs. If a measurement tells us the mass and radius of a NS, we can equally con- 
clude that in GR, some EOSs are consistent with this measurement, or for a specific 
EOS, some gravity theories can pass the test. Thus, in order to constrain gravity 
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theories, one needs to consider various reasonable EOSs, or one has to break the 
degeneracy between gravity theories and nuclear physics. 

One way to break such a degeneracy is to find relations that are not sensitive to 
the EOS but still depend on the gravity theory. In GR, such relations do exist. One 
famous example of such relations is the “J-Love-Q” relation between the dimension- 
less moment of inertia, the tidal Love number and the quadrupole moment of NSs 
(see, e.g. Fig. 1 in Ref. [125]). The deviations from the common relation caused by 
variations in EOS can be smaller than 1%. As the tidal Love number can be mea- 
sured by GWs from colliding NSs as mentioned before, measuring the moment of 
inertia of NSs thus becomes very important. 

Many efforts have been done to understand why there are such universal relations 
existing in GR and many other gravity theories [126]. The possible physical picture 
is some approximate symmetries emergent when the stars become more and more 
compact and they are responsible for the universality [126]. For example, the iso- 
density self-similarity leads to a constant eccentricity profile which will minimize 
the energy of the system [62]. There are also some mathematical approaches show- 
ing that the “I-Love” relation is perturbatively insensitive with respect to changes in 
the polytropic index of the EOS around the incompressible limit [18]. 

There are many other universal relations known today, for example the relation 
between the Love number and the f-mode frequency, the no-hair-like relations, the 
I-C relation and so on [126, 44]. Testing strong-field gravity with universal relations 
seems promising, though there are still some subtle issues that need to be studied fur- 
ther [102]. In principle, to test gravity with universal relations, one needs to measure 
several quantities independently for the same NS, for example, the Love number 
and the moment of inertial. But usually, different quantities are proper to measure 
in different systems, which correspond to different NSs with different masses. One 
way to go through this problem is to measure those quantities for many different 
NSs with similar masses [102], for example around 1.4Mọo. Then one can combine 
these results to a universal relation by doing an interpolation on the mass. Another 
possible method is to establish universal relations that relate quantities in different 
masses [92]. 


1.5 Tests of Alternative Gravity Theories 


Though Einstein’s GR has passed all present tests with flying colors, there is still 
great enthusiasm for modified gravity theories [1 18, 10, 117, 100]. Those alternative 
gravity theories may arise as low-energy effective theories of some quantum field 
theories or motivate from dark matter and dark energy problems. As they provide 
the explicit form of how can a theory deviate from GR, studying these modified 
gravity theories is also important and intuitive in testing gravity. In the following, 
we will describe some basic aspects of modified gravity theories by illustrating two 
concrete examples: scalar-tensor gravity and massive gravity. 
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1.5.1 Scalar-tensor theory 


Scalar-tensor gravity theories are the simplest and mathematically well-posed ex- 
tension to GR by including some scalar fields that are non-minimally coupled in the 
Lagarangian of gravity. A general class of scalar-tensor gravity theories with one 
extra scalar field ® have the following action in the physical frame (or the so-called 
Jordan frame) [25] 


g dx z 
_ - B a E 
= CA c V—8[F(®)R-Z(®) gs" dy Bay -U (D) 
+Sm[Wins Suv] A (1.40) 


where G, denotes a bare gravitational coupling constant and Sm[Wm; uv] is the ac- 
tion of standard matters; R and g are the Ricci scalar and the metric determinant 
calculated from the “physical metric” guy, respectively. In this kind of theories the 
WEP is satisfied because all the matter variables Ym are coupled to the same metric 
uv. It is often convenient to rewrite this action in a canonical form by redefining ® 
and Syy via [22] 


Suv = F(P) Suv, (1.41) 


3F2(@) 1Z(®))'” 
(= + [ae $ F2? (®) a 2 F(®) 


(1.42) 


It leads to an action in the Einstein frame 


paat Ue ara 28!” ð 90,9 —V(@)] 
dened c V Ee ee SEE YP 
+Sm[WinsA7(Q) Suv] - (1.43) 


Now R, and g, are calculated from the new metric guv» and the potential term 
V (p) = F~*(®)U(®). (1.44) 


In this form the matter fields are coupled to the conformal metric A*(@) guv instead 
of giy itself, where the conformal coupling function 


Alo) =F (®). (1.45) 


As the matter fields couple directly to the metric uy, this metric is the one that can 
be measured by laboratory clocks and rods and it is called the physical metric. While 
the canonical representation is often more convenient for calculation as it decouples 
the spin-O excitation of @ and the spin-2 excitation of Euv: 

For different choices of the coupling function A(@) and the potential V(@), one 
gets different classes of scalar-tensor theories. Note that one can always work in 
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some appropriate units that the asymptotic value of A(@) = 1, which means in these 
units, the Einstein metric and the physical metric asymptotically coincide. 

The field equation can be derived from Eq. (1.43) by taking variations with re- 
spect to giy and @ [123] 


1 1 
Ruy = 29u Pv@ + 82G, (ti = Tew) + 58uvV (9), (1.46) 
1 dV(o) 
= —AnG, pao 1.47 
o nGa (P)T; + A dp (1.47) 


where T” = 2(—g4) 1 ES m/ 88v and a (p) = dInA(@)/d@. The quantity a(@) 
plays the role of the coupling strength between the scalar field and ordinary matter, 
and in general, it can be field-dependent. 

Since any function a(@) =InA(@) can be expanded when the scalar field is weak 
as 


1 
a(p) = a0(P — Po) + Bol — Po)? ++, (1.48) 


it has been shown that, for V(@) = 0, in the appropriate units that a(o) = 0, all 
measurable quantities depend only on the first two coefficients œo and Bp of the 
expansion at the 1 PN level for weakly self-gravitating objects [25]. In the PPN 
formalism [118], only two parameters y?PN and BPPN deviate from GR (in which 
yPPN = BPPN — 1), and they are given by 


202 a Bo 
PN 0 PPN 0 


One can also find that the effective gravitational constant between two bodies is 
G=G,(1+ a) in the weak field. 

A famous theory, which also may be regarded as the simplest case for scalar- 
tensor gravity, is the Brans-Dicke theory [14], where there is a varying gravitational 
constant. In the scalar-tensor formalism, the Brans-Dicke theory is equivalent to 


a(@) = Q@-+const., V(g) =0, (1.50) 


and in this theory there is a field-independent coupling strength, a@(@) = a. Solar 
system experiments including light deflection, time delay and frequency shift mea- 
surements made with the Cassini spacecraft have set limit on the coupling strength 
of the scalar field that Qe < 1 x 1075 [88, 91, 11], which tightly constrains the de- 
viations from GR in the weak-field region. 

Things can be much different for systems including strongly self-gravitating 
objects like NSs. It was found that a wide class of scalar-tensor theories can ex- 
hibit nonperturbative behaviors called spontaneous scalarization for strong-field ob- 
jects [26, 34]. As those theories can pass all the weak-field gravity tests while pro- 
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viding large deviations from GR in the strong-field region and still keep good math- 
ematical behaviors, they largely enhance the interests in studying such theories. 

For strongly self-gravitating objects, one can define the effective scalar charge 
of an object A by œa = d lnm4 /9 Q0, or equivalently by the asymptotic behavior of 
the scalar field [25] @ = Po —G.maaa/r+O (1/r?). The Newtonian gravitational 
coupling between two bodies A and B can be described by the effective coupling 
constant, Gag = G..(1+ 04a). 

A simple theory that can have spontaneous scalarization is called the Damour- 
Esposito-Farése (DEF) theory, which is widely studied [26, 22]. It can be described 
by A(@) = eP®/2 and V(@) =0, and the asymptotic value @ of @ here we set to 
be zero. With arbitrary choice of the asymptotic value ọọ, the DEF theory is just 
a simple choice with œo = P o, Bo = P in Eq. (1.48). In this theory, the coupling 
strength a&(@) = BQ is a linear function of 9. 

The name spontaneous scalarization used by Damour and Esposito-Farése comes 
from the analogy to the phenomenon of spontaneous magnetization of a ferromag- 
net [27, 94]. In the DEF theory, the trivial case that @ = 0 is always a solution to 
the field equations. In fact this solution is identical with that in GR. In the weak- 
field region, as we have set @o = 0, the DEF theory reduces to GR so it can pass 
any weak-field gravity tests that GR can pass. But for a strongly self-gravitating 
body, a large scalar charge may be produced spontaneously, like a ferromagnet will 
have spontaneous magnetization at low temperature. The solution with a non-vanish 
scalar charge is energetically favored when the star’s baryon mass exceeds the criti- 
cal mass [26, 34]. 

One can also understand this phenomenon by analyzing the instability of the 
perturbation in the scalar field [86]. To the linear order in @, the field equations 
become 


Riv = (1.51) 
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ọ = —4aT Bo. (1.52) 


For matter described by a perfect fluid, we have T = 6 +3), which is negative for 
ordinary matter. Thus the right-hand side of Eq. (1.54) is negative for a negative p, 
and the scalar field can suffer the tachyonic instability for a sufficiently negative ß. 
Interestingly, it has been shown in some works that for highly compact NSs with 
certain EOSs, T can be positive, which leads to a spontaneous scalarization for 
p > 0 [80]. 

For those scalar-tensor theories with spontaneous scalarization, the weak-field 
tests cannot constrain their parameter space efficiently, but they still can be con- 
strained from strong-field tests. In scalar-tensor theories, the energy loss of a system 
can be carried by both the spin-2 GWs and spin-0 scalar waves to infinity. For those 
theories with spontaneous scalarization, NSs can carry a large number of scalar 
charges, which will cause a much enhanced gravitational dipolar radiation for a 
scalarized NS in a binary. To the leading order, it will contribute an additional decay 
rate of the binary orbital period via [27] 
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Fig. 1.7: The dashed black curves show the effective scalar charge in the DEF theory 
with |æo| = 1075 and By = —4.8, —4.6, —4.4, —4.2 from top to bottom repectively. 
The EOS is assumed to be AP4. The observational bounds from different binary 
pulsars and a GW event are denoted by the triangles [97]. 


přpole = i 
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20Gx 5/920 mam 
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For some binary pulsars, while a NS has a significant scalar charge, its companion, 
for example a WD, can only take a vanishingly small scalar charge ag — 0 as it 
is a weak-field object. Therefore, NS-WD binaries can be the most sensitive probe 
to constrain the parameter space of spontaneous scalarization [42]. There is also 
study showing that for NS-NS binaries with significant difference in the masses of 
the components can also be excellent laboratories [130]. Besides that, it has been 
proven that the no-hair theorem still holds in many scalar-tensor theories, which 
means that the BHs in those theories cannot carry a scalar charge. So that although 
the binary BHs in those theories will not excite dipole radiation, the NS-BH binaries 
can still provide abundant information about the underlying gravity theory [69, 98]. 
Figure 1.7 illustrates the constraints on the effective scalar charges from the grav- 
itational dipole radiation of seven binary pulsar systems and one GW event [97, 
130]. Note that in this figure, a specific EOS, AP4, is used to draw the curves. In 
principle, the uncertainty in the NS EOS will entangle with the gravity test. Nev- 
ertheless, with enough well-measured binary pulsar systems, one can set a limit on 
the whole mass range for NSs, and studies have shown that for each reasonable 
EOS, the parameter space for spontaneous scalarization in the DEF theory is now 
very small. In Refs. [99, 131, 46, 130], by using the Bayesian parameter estimation, 
the possibility of the spontaneous scalarization for an effective coupling larger than 
107? for DEF theory is basically closed, no matter what the underlying EOS of NS 
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As the parameter space for spontaneous scalarization in the initial DEF theory 
is largely constrained by pulsar timing observation, a natural way to avoid the con- 
straints while keeping the phenomenon of spontaneous scalarization is to consider 
massive scalar-tensor theories [86, 34]. For example one can consider the DEF the- 
ory but with V (p) = 2m? ọ?, where m is the mass of the scalar field. The linearized 
field equation of the scalar field now becomes 


ọ = (mP —4nTB)oQ. (1.54) 


It shows that the mass term of the scalar field will suppress the existence of sponta- 
neous scalarization [86]. Nevertheless, the spontaneous scalarization can still hap- 
pen for massive scalar-tensor theories. But with a mass term, the asymptotic value 
of the scalar field is suppressed exponentially in a Yukawa fashion, which means by 
definition the scalar charge of a scalarized NS vanishes. Thus no significant dipole 
radiation will go to infinity. What is more, when the Compton wavelength of the 
scalar field is smaller than the orbital separation of the binary, the modification with 
respect to GR in the orbital dynamics will also be suppressed exponentially, which 
causes the pulsar timing observation hardly setting strong constraints in some pa- 
rameter space of massive scalar-tensor theories. Fortunately, this kind of theory can 
still be constrained via the tidal deformability measurement from GWs (see e.g. 
Refs. [123, 55]). This is true for the reason that a scalarized NS can have large 
difference in radius compared to that in GR, while the tidal deformability is very 
sensitive to the star’s radius. In this sense, combining pulsar timing observation and 
GW observation can probe a large parameter space of scalar-tensor theories [99]. 

In recent years, some variants of scalar-tensor theories triggered great enthusiasm 
as they can provide solutions for scalarized BHs, in contrast to the no-hair theorem 
(see Ref. [34] for a review). One example of those theories is the scalar-Gauss- 
Bonnet theory, which includes a topological Gauss-Bonnet term, 


G = RopysR*? —4Rap RP +R, (1.55) 


which is coupled to the scalar field [35, 104, 124]. Scalarization of NSs in such kinds 
of theories can be very different from those in the DEF theories. It also provides 
a new field where observations of compact objects, including NSs and BHs, are 
crucial to revealing the strong-field information of gravity. More examples of other 
variants of scalar-tensor theories that pass current binary-pulsar tests can be found 
in Ref. [34]. 

As already discussed, the theories of gravity can sometimes degenerate with the 
EOS of NSs. One way to break the degeneracy is to find those universal relations 
that do not depend on the EOS. We have shown that such universal relations do exist 
in GR, and in fact, this can still be true for some modified gravity. For some theo- 
ries without spontaneous scalarization, universal relations in GR, like the I-Love-Q 
relation, are still valid. That means in modified theories, these relations also only 
weakly depend on the EOS of NSs [126], but they can be very different from their 
GR counterpart. Things can be more interesting when spontaneous scalarization 
kicks in. The universal relations in GR can be EOS dependent for those theories 
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with spontaneous scalarization, because in general, the starting point of the spon- 
taneous scalarization is different for different EOSs. Once a NS is spontaneously 
scalarized, the relations can differ from their GR counterparts largely. Thus for the- 
ories with spontaneous scalarization, we can see that the relations for different EOSs 
branching from the GR one from different starting points. For some theories, like 
the DEF theory, those different relations for different EOSs will tend to form a new 
universal relation different from GR when the scalarization becomes stronger. How- 
ever for some theories like the massive Gauss-Bonnet theory, this may or may not 
happen depending on the theoretical parameters [124]. 


1.5.2 Massive gravity theory 


From a modern particle physics perspective, the related quantum particle for gravity 
in GR is a massless spin-2 particle, the so-called graviton. However, massive gravity 
theories, which maintain the notion that gravity is propagated by a spin-2 particle 
but consider this particle to be massive, have aroused wide interest. Such a massive 
gravity theory may be used to explain the acceleration of the Universe and dark 
matter [52]. Probing the upper limit of the graviton mass is a basic topic in physics. 
Constraints from the pulsar timing technique may come from two aspects: one is the 
radiative tests from binary pulsars, and the other is the orbital motion of a pulsar- 
supermassive BH system. 

Unlike the dispersion relation tests from the LIGO/Virgo/KAGRA GW events [2], 
the constraints on the graviton mass are generally theory dependent in pulsar tim- 
ing, therefore we shall introduce some specific massive gravity theories. One of the 
examples is the cubic Galileon theory, which is a simple model with the Vainshtein 
screening mechanism [31]. The action including matter is 


1 hEVT, 1 T 
S= Ja uncon: at (gn) (1+ z.) + g | ; 


2Mp, 4 3A3 2Mp\ 
(1.56) 


where huv = Suv — Nuv, (Eh)uv = -4 huv +++- and T is the trace of the stress- 
energy tensor. In the action, Mp; is the Planck mass, and the strong coupling scale of 
the Galileon sector A is related to the graviton mass mg via A? = m;Mpi. The field 
equation for ms and huy are [31] 


1 1 

——Tyy = — -0h 1.57 
Mp, uv 2 UV; ( ) 
l 7 9,|-2a"s, (1+ 0r + Lan (an)? (1.58) 
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For a static point source with mass M, the solution for 2 can be written as Vr, (r) = 
fE (r), where 
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and the Vainshtein radius associated with the mass M is rą = (M /16A 3Mp1) 


1/3 
(m / 16m2M? ) . Within the radius r,, the fifth force from the scalar field is 


strongly suppressed, and the theory reduces to the canonical gravity. With such a 
screening mechanism, this theory can avoid the stringent constraints from the Solar 
system while it still can deviate from GR at a large scale, thus providing changes to 
the cosmological evolution. 

Though the screening mechanism tends to suppress the deviation from GR at 
high-density regions, the cubic Galileon theory still predicts difference in gravita- 
tional radiation [31]. For a binary system with a characteristic velocity v, the grav- 
itational radiation is less suppressed by a factor of ~ y3/2 compared with its fifth- 
force counterpart. As shown by de Rham et al. [31], the extra radiative channels in 
these theories include monopolar radiation, dipolar radiation and quadrupolar radi- 
ation, and for different systems, for example binary pulsar systems with different 
orbital periods and orbital eccentricities, the dominate radiation channel can be dif- 
ferent [101]. Recent study with 14 well-timed binary pulsar systems has set a tight 
constraint on the graviton mass in the cubic Galileon theory, mg < 2 x 10-8 eV /c?, 
at 95% confidence level [101]. 

Another example we discuss here is the so-called Yukawa gravity, which is a 
widely used phenomenological form that contains the possible deviations from the 
Newtonian gravity via a modified potential, 


pe ae Fr: s 
Ees T (1+ae i (1.60) 
Part of the Newtonian potential is suppressed exponentially, with a strength of a 
and a length scale of Ay. For œ = 0 or Ay = œ, the Yukawa potential will reduce 
to the Newtonian one, and the graviton mass here can be expressed as mg =h/Ay. 
Such a Yukawa gravity can naturally arise in the Newtonian limit of f(R) gravity 
with a general action [17] 


s= [axv=8 |ræ+ 2, | (1.61) 
C 


Assuming that the function f(R) is Taylor expandable and take the O(2) approxi- 
mation, f(R) = fo + fiR+f2R? +--+ with constants f; (i = 1,2,3), the gravitational 
potential defined via gy = 1+ 2® will take the form of Eq. (1.60), and the Newto- 
nian potential is obtained only when f(R) = fiR. 

Tests of the Yukawa gravity have been carried out at a wide range of length 
scales with different systems and methods, including terrestrial laboratory exper- 
iments [83], LAser GEOdynamic Satellite (LAGEOS) [85], Lunar Laser Ranging 
(LLR) [119], Solar system planetary orbits [107], S2 star orbit [129] and observa- 
tions of elliptical galaxies [82]. The reason for taking experiments at such a wide 
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Fig. 1.8: Constraints on the strength of the Yukawa interaction œ as a function of 
the graviton mass mg from different experiments (figure courtesy of Y. Dong) [36]. 
Shaded regions show the excluded parameter space at the 95% confidence level. 


coverage is that for the Yukawa gravity theory, only when the system scale L is 
similar to the length scale Ay it shows the largest deviation from Newtonian grav- 
ity. Thus timing a pulsar orbiting around the supermassive BH in our galaxy center 
may provide a unique opportunity to test the Yukawa gravity theory at a new length 
scale [36]. For a pulsar orbit with a semi-major axis a and an orbital eccentricity e, 
it will experience an additional periastron advance from the Yukawa potential at the 
rate of [128] 


1 maV1—e2 a 
@ = i (1.62) 
P, l+a Ay 


to the first order in a/Ay. If the pulsar orbit has a semi-major axis a ~ Ay, periodic 
effects such as the deformation of the orbit may also be detectable by pulsar timing. 

Figure 1.8 shows the projected constraints from a S2-like like pulsar, on the 
strength @ at 95% confidence level as a function of mg [36], from where one can 
see that the pulsar-supermassive BH system can contribute a unique test of this the- 
ory that is not covered by other experiments. For a given strength of the Yukawa 
interaction œ one also can obtain the related constraint on the graviton mass, and 
from Eq. (1.62) one can see that for œ changing from | to ~, the periastron advance 
rate only changes by a factor of 2, thus the constraints on mg will not change too 
much. Simulation results by Dong et al. [36] showed that the expected constraints 
on graviton mass can reach mg < 10-74 eV /c? for a 20-yr observation of a S2-like 
pulsar around Sagittarius A*. 
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1.6 Summary 


In this chapter, we present a pedagogical introduction to some basic aspects of pulsar 
timing and using it in testing gravity theories. As pulsar timing provides very high 
accurate measurements of a strongly self-gravitating object in orbit, it can put strict 
constraints in the strong-field region of gravity. Moreover, the observational preci- 
sion of the pulsar-timing parameters generally improves with the observational time 
span Tops. For example, the precision in the orbital decay parameter, P,, improves 


as Ta? ? Thus long-term observations of interesting systems like the Hulse-Taylor 
pulsar and the double pulsar can provide us with new results from time to time until 
the precision is dominated by other effects (e.g. the Galactic potential model). Fur- 
thermore, the sensitivity of radio telescope is improving at the same time, so the real 
improvement of timing parameters is faster than the theoretical power law predic- 
tions. New telescopes like the FAST telescope in China [74] and the next-generation 
radio telescope, the Square Kilometre Array in South Africa and Australia [116] can 
provide plenty of new results. New systems like pulsar-BH binaries or pulsars near 
the Galactic center may be discovered, and they may lead to new tests of gravity and 
give better results [71, 69, 98]. 

As we mentioned in the chapter, combining pulsar timing observation with other 
observations can probe a large theory space of both gravity theories and nuclear 
physics. GW observations [2] and BH shadows [3] all provide unique information of 
gravity in different gravity regimes. For GWs, with the third generation of detectors 
that are under construction, almost all binary NS mergers in the Universe can be 
detected, and the BH shadows and X-ray tests provide the possibility of measuring 
the BH metric in the near zone. Combining with proper universal relations, one 
can break the degeneracy between gravity theories and nuclear physics, thus giving 
EOS-unbiased tests of gravity. 
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